Repulsively induced photon super-bunching in driven resonator arrays 
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We analyze the non- equilibrium behaviour of driven nonlinear photonic resonator arrays under the 
selective excitation of specific photonic many-body modes. Targeting the unit-filled ground state, we 
find a counter-intuitive 'super bunching ' in the emitted photon statistics in spite of relatively strong 
onsite repulsive interaction. We consider resonator arrays with Kerr nonlinearities described by the 
Bose Hubbard model, but also show that an analogous effect is observable in near-future experiments 
coupling resonators to two-level systems as described by the Jaynes Cummings Hubbard Hamiltonian. 
For the experimentally accessible case of a pair of coupled resonators forming a photonic molecule, 
we provide an analytical explanation for the nature of the effect. 



Introduction - Coupled resonator arrays (CRAs) offer 
the intriguing possibility of realising strongly-correlated 
many-body quantum states of light. Early work on CRAs 
assumed idealised, lossless arrays, and focussed in partic- 
ular on equilibrium quantum phase transitions in these 
structures. However, near-future photonic devices will 
necessarily operate under driven-dissipative conditions 
on account of unavoidable photon loss, thereby serving 
as natural platforms for the exploration of novel non- 
equilibrium many-body photonic effects [IHS]- Our un- 
derstanding of these systems is in its infancy, making 
it desirable to concretely connect the non-equilibrium 
properties of CRAs with their more familiar equilibrium 
structure. To this end, there have been recent efforts 
to identify signatures of the equilibrium quantum phase 
transition as originally proposed in p^OHTS] which survive 
under lossy dynamics. 

We propose an alternative scheme to chart dif- 
ferent regions of parameter space and connect non- 
equilibrium observables to the underlying Hamiltonian 
properties. We envisage a resonator array driven to a 
non-equilibrium steady state (NESS) by external lasers, 
with the laser frequency chosen such that the unit-filled 
equilibrium ground state with on average one particle 
per site is selectively addressed and populated. Fea- 
tures arising from the details of the non-equilibrium op- 
eration appear in collected emission statistics, including 
a counter-intuitive many-body repulsion-induced bunch- 
ing of the emitted photons, the magnitude of which is 
controllable via tuning Hamiltonian parameters. Novel 
super bunched light sources far exceeding the bunching 
of thermal photons may find important applications in 
ghost imaging technologies [19 and all-optical simulation 
of two-photon correlations in quantum walks [20j. 

In this work we first focus on a minimal-sized two-site 
resonator system. Such dimers or a photonic molecules 
are expected to be experimentally viable in the near fu- 
ture in different technologies ranging from semiconduc- 
tors to Circuit QED. We initially study the system for 
resonator nonlinearities of the repulsive Kerr-type, as 
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FIG. 1. (a) Schematic of the driven-dissipative Bose Hubbard 
model, featuring local coherent driving, photon tunnelling 
and a purely photonic Kerr nonlinearity. (b) The driven- 
dissipative Jaynes Cummings model with effective photon 
nonlinearity generated by couplings to two-level systems, (c) 
Diagram showing the bare basis (solid black) and eigen- 
frequencies (dash dotted red) of the driven system for the 
minimal system of M = 2 resonators. The laser is tuned so 
that two laser photons (vertical red arrows) are capable of 
promoting the system to the lowest-lying 2-photon state. 



shown in Fig. [T] (a), to illustrate our driving scheme. 
We then analyze the Jaynes- Cumming type encountered 
when resonator modes interact with embedded effective 
two- level systems [21 , as in Fig.[l](b). We note here that 
we have previously investigated the validity of modelling 
the JCH with the simpler, single species Bose-Hubbard 
(BH) model [22 , and accordingly also demonstrate that 
the super-bunching signature persists under a JCH de- 
scription. Moving beyond this minimal 'array', bunched 
emission is also demonstrated in near future experimen- 
tally accessible mesoscopic-sized systems [23, • 

System - We consider a one-dimensional chain of M 
coupled single-mode optical resonators under periodic 
boundary conditions. Each resonator of frequency is 
coherently coupled to its two nearest-neighbors. Exter- 
nal lasers coherently drive each resonator in-phase with 
amplitude 1]. In a frame rotating at the laser frequency 
uji,^ the system Hamiltonian is: 
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Here is the Kerr nonlinear strength, J is the photon 
hopping rate, Vt is the photon driving strength, and (j, j') 
denotes nearest neighbour resonators. The operators CLj 
are the photon destruction operators for the photon mode 
in resonator j. The detuning of the driving laser fre- 
quency from the bare cavity frequency is Ac = co'c — cj/,. 

Markovian photon loss processes from each cavity are 
incorporated via a quantum master equation formalism 
for the evolution of the system density matrix p, p = 
C^^\p], where 

M 

The dissipative part of the dynamics is described by 
V^[p] =OpOt-i(OtOp+pOtO). The NESS of the sys- 
tem is described by the density matrix which satisfies 
J^[pss] = and observables are measured with respect to 
this state, (0)ss = Tr(Opss)- 

Two-resonator ^dimer^ - We begin by analyzing the 
simplest possible driven resonator 'array' consisting of 
just M = 2 resonators, which serves to illustrate clearly 
our scheme for accessing the unit-filled ground state. 
Figure [l] (b) shows the low-lying eigenstructure of the 
undriven Hamiltonian of Eq. ([T]) for M = 2, and our 
driving scheme. The two one-photon eigen-frequencies 
are the symmetric (+) and anti-symmetric (-) Bloch 
modes, and lie (in the bare frame, with Q = 0) at 

= ujc T with corresponding eigenstates 
The two-photon eigen-frequencies are uJq ^ = 2uOc -\- 2U^ 
uj^^^ = 2uJc ^ U T Vt/^ +4J2, with eigenstates |2o), 
|2±), respectively. The unit-filled ground state is of fre- 
quency uj^ . This mode undergoes a qualitative change 
between the extreme limits of a localized state, charac- 
terized by vanishing on-site photon number fiuctutation 
Var(a^aj) ^ for /7 ^ J, to a coherent superposition 

state with Var (a^a^) ^ ^ for <C J. For increasing sys- 
tem size, the behavioural transition of the ground state 
becomes sharper, approaching the celebrated Bose Hub- 
bard Mott-insulating to superfiuid phase transition in the 
infinite system limit [25 . 

To selectively populate the unit-filled ground state of 
Eq. ([1]), we set the driving laser frequency such that two 
laser photons are resonant with the lowest-lying two pho- 
ton mode, i.e. 2ujl = cj^ , implying a laser detuning 




FIG. 2. (a) Emitted photon statistics as a function of J and 
JJ for a dimer (M = 2) of resonators driven according to 
Eq. ([3]). Black curve: analytic result for g^^^ — 1. Black dot: 
the point ( Jc/7p, t^c/7p) where bunching sets in. (b) 'Slices' 
along the dotted lines in (a), above and below the critical 
hopping amplitude at which the bunching feature appears, for 
decreasing driving strengths, down to the infinitesimal limit 
(solid black lines). 

Fixing Ac = Ac(J, f/) for a given hopping and nonlin- 
earity, we now examine which features of the underlying 
Hamiltonian mode structure leave fingerprints on exper- 
imentally accessible photonics observables. In our non- 
equilibrium setting the photon number is not an integer 
and the photon number variance is not an informative or- 
der parameter. To compensate for these non-equilibrium 
effects, we instead focus on the local zero-time photon 
correlation function g^^^ = g^^^ = {a^-a^-ajCij) / {a^-CLj)'^ ^ 
a statistical quantity directly accessible in CRA setups 
through standard methods like homodyne detection. We 
note that g^'^^ measurements may be particularly valuable 
in weakly-driven systems, where the excitation number 
may be very small, but normalised statistical quantities 
may be collected over longer times. 

Figure [2] (a) shows g^'^^ measured in the NESS of a 
BH dimer pumped at the laser detuning of Eq. ([3|, for a 
range of tunnelling rates and nonlinearities (J, /7). The 
diagram is broadly divided into three regions, defined 
by Poissonian {g^'^^ ^ 1), anti-bunched {g^'^^ < 1) and 
bunched {g^'^^ > 1) statistics. Notably, there is a critical 
coupling rate Jcrit between the resonators below which 
bunching does not occur for any value of nonlinearity, 
suggesting that the bunching arises from a co-operative 
many body effect in the NESS. Figure [2] (b) shows the 
qualitative difference in the behaviour of the correlation 
function above and below this critical point. For compar- 
ison, an isolated resonator driven at it's single-particle 
(unit-filled) resonance never exhibits bunched signatures 
(see dash-dotted line in lower panel of Fig. [2] (b)). 

At low nonlinearities U <C 7p, the dimer is driven 
at the frequency of the zero-momentum Bloch mode, 
and its response is approximately linear. The NESS 
is a coherent state, inheriting Poissonian statistics with 
^^^^ = 1, and average population {ojCtj) = (2^7/7^)^, for 
all J. At the other extreme, taking the hardcore limit 
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U ^ oo, no more than a single photon per resonator 
can be injected, and the problem is reduced to two cou- 
pled two-level systems whose emission is completely anti 
bunched (^^^^ = 0), with mean excitation number per 
resonator limf/^oo(^^^) = x{2x + l)/((2x + 1)^ + xy), 
where x = {2ft/ jp)'^ and y = [J/Q)'^. 

Away from these extreme limits, in the region of pa- 
rameter space where J and U are comparable, the emit- 
ted light is bunched, for hopping rates larger than a criti- 
cal rate J > This is counter-intuitive, as we are prob- 
ing a two-photon state with significant repulsion favour- 
ing separation, and yet we find an enhanced probability 
of photons being emitted together, relative to the statis- 
tics of the driving. 

We derive features of this bunched region by consider- 
ing the limit of infinitesimal driving. Figure |2] (b) shows 
that the correlation function approaches a limiting be- 
haviour as fl/^p 0, an observation confirmed by per- 
turbatively expanding the elements of the NESS density 
matrix p^s in increasing powers of the driving strength, 
then solving for g'^'^^ [8 . The equations thus obtained are 
physically opaque, however further progress can be made 
by making a pure-state ansatz pss = |^ss)(^ss|7 valid in 
the low driving regime [7]. The state |^ss) is found as the 
stationary state (corresponding to the zero eigenvalue) 
of the effective Hamiltonian Ti^^^ obtained by replacing 
Ac Ac — ijp/2 in Eq. ([T]), which may be interpreted as 
the Hamiltonian governing a single quantum trajectory, 
with a vanishing probability of a quantum jump ensured 
by taking the limit ft/jp 0. 

Considering only the lowest excitations and exploit- 
ing the symmetry of the two-site system, we set |^ss) = 
CoolOO) + Ci(|01) + |10)) + Ciilll) + C2(|02) + |20)), 
from which the correlation function follows as ^^^^ = 
|C'2p/|Ci|^. We obtain the minimum resonator coupling 
for which bunched statistics appear as Jc = J dip — 

^(3 + 2>/2)/4, and the nonlinearity at this point IJc = 

Uc/lp = \/7c as marked in Fig. [2] (a) (see supplementary 
material for further details). The transition from super- 
to sub-Poissonian statistics (i.e. where g^'^^ = 1) is found 

to occur at J ~ -^/7/2 for J ^ Jc^ while for very large 
but finite /7/J, the exact solution may be simplified to 




Figure [3] offers physical insight into this phenomenon, 
showing the emission spectrum of the system as calcu- 
lated from the Fourier transform S{uj — ul) of the on- 
site steady-state auto-correlation function S{t) = {a^{t-\- 
r)d{t)), as a function of increasing nonlinearity. The res- 
onator coupling is sufficiently large {J/jp = 10 > Jc) to 
observe bunched emission (top panel of Fig. [3] (a)). At all 
nonlinearities, the spectrum is dominated by two bright 
features. These correspond to decays from the lower and 
upper one-particle states to the vacuum |0), labelled 
lines A and B respectively. 



Weaker features are also present, which do not signif- 
icantly affect the steady state photon populations, but 
may strongly modify statistical quantities such as 
Line C in Figs, [s] (a) and (b) corresponds to emission 
from the highest two-particle state |2+) to an intermedi- 
ate level, as drawn in Fig. [3] (c). In the vicinity of the 
crossing of emission lines B and Cat[//J=(9 — vT7)/4, 
the population in |2_) reaches a maximum, as photons 
emitted as part of the line B decay process may transfer 
population instead to |2_). Further, calculating projec- 
tions into the bare basis, we find that there is a greater 
probability of finding the two photons of the mode |2_) 
in the same resonator than distributed between them, 
relative to the driven |2+) mode. This is reflected in 
the enhanced probability of simultaneous emission of two 
photons (^^^^ > 1) around this crossing. In contrast the 
mode 1 2+) favours delocalising its two photons relative 
to the statistics of the driving laser. Thus, we observe 
either approximately coherent, or anti-bunched light in 
all regions of parameter space except in the vicinity of 
the crossing of lines B and C. For resonator couplings J 
less than the line width 7^ of the resonators, the global 
physics of the system resembles that of an isolated nonlin- 
ear resonator driven at it's resonance frequency, such that 
anti-bunching is always expected - in spectral terms, the 
crossing of lines B and C is hidden inside the coalesced 
lines A and B. 




FIG. 3. (a) Lower: Emission spectral function relative to the 
laser frequency \S{uj — ul)]'^ with increasing nonlinearity U at 
a fixed hopping J/jp = 10. At our weak driving {ft/^p = 0.3), 
the spectrum is dominated by transitions between the ground 
and one-particle manifold (bright lines A and B). Weaker fea- 
tures involving the two-particle manifold are also present, for 
instance the dashed curve highlighting line C. Upper: the 
zero-time correlation function g^'^\ (b) Magnification of the 
region inside the white square in (a), calculated at the higher 
driving fl/'^p — 1 to highlight the crossing of the features 
labelled Lines B and C. (c) Transitions involved in spectral 
lines B and C, with only the relevant modes drawn. 

Larger systems - We now investigate how the corre- 
lations presented in Fig. |2] (b) evolve as the system size 
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increases, continuing to drive the commensurately- filled 
ground state. An analytic approach valid for arbitrary 
Hamiltonian couplings beyond M = 2 resonators is in- 
tractable. Instead, we numerically calculate the eigen- 
vector of the effective Hamiltonian l-L^^^ with eigenvalue 
closest to zero (taking a series of successively weaker driv- 
ings Vt/^pio ensure convergence of observables). Exploit- 
ing the translational invariance of systems with periodic 
boundary conditions allows us to access systems of up 
to M = 7 resonators while retaining three photons per 
resonator in simulations. Rigorous quantum trajectory 
calculations based on the matrix product state repre- 
sentation and the time evolving block decimation algo- 
rithm [26H28] performed at a finite driving strength (see 
Supplementary material) broadly agree with the results 
obtained via numerically exact diagonalisation of i-L^^^ • 
However, they do indicate that the precise features of the 
bunching region do depend on drive strength, as already 
observed for M = 2 in Fig. |2] (b). 

Figure [4] shows the evolution of the counterintutive 
bunched region for increasing system sizes up to M = 7 
resonators. We see a reduction in the magnitude and 
range of interaction strengths for which bunched light 
is observed as the system size increases. The bunching 
region is seen to retreat up the J- axis, while smaller in- 
teraction strengths U are necessary to induce the bunch- 
ing. This explains the reduction in the magnitude of the 
effect observed for cross-sections at constant resonator 
coupling, as in Figure |4] (a). 
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FIG. 4. (a) Correlation function evaluated as a function 
of increasing nonlinearity at fixed resonator coupling J/^p = 
10^ Note the M = 2 results are not included as the bunching 
in this case is significantly larger, (b) Extent of the bunched 
region in (J, U) parameter space, as measured from U — [/lhs 
corresponding to the peak g^'^\ to U — ?7rhs at which the 
correlations change from bunched to anti-bunched, for a range 
of resonator couplings J. 



Photon statistics in Jaynes Cummings arrays - Near 
future circuit QED systems will most probably re- 
alise a few-photon resonator nonlinearity via a Jaynes- 
Cummings interaction with embedded two-level systems 
[lOj . as described by the Jaynes Cummings Hubbard 
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{^cd]dj + (A, - A)(7+a^. 
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Here A = — denotes the difference between the 
resonators' frequency and the TLS transition frequency, 
and (j^ denote TLS raising and lowering operators. The 
JCH is known to possess a localized-delocalized transi- 
tion as either the hopping J is increased, or the Jaynes- 
Cummings parameter A is made more negative. This 
transition is similar in some respects to the phase transi- 
tion of the BH model, though also differs in fundamental 
ways on account of the different nature of the systems's 
intrinsic excitations (bosons and polaritons, respectively) 

[HHTHIES]. 
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FIG. 5. Steady state observables for a two-site Jaynes- 
Cummings array, driven at its lowest two-particle resonance. 
Parameters: ^/7p = 10. (a) Zero-time photon correlation as a 
function of resonator coupling J and atom-resonator detuning 
A, playing the part of an effective photon nonlinearity. (b) 
Spectral function |5'(c(;)|2 evaluated along the dashed line in 
(a), again showing a crossing of emission lines promoting pho- 
tons to a state with enhanced probability of bunched emission 
(red arrow is the driving laser frequency). 

Figure |5] presents evidence that the mechanism under- 
lying the bunched emission discussed above in the con- 
text of a the driven Bose- Hubbard model persists in this 
qualitatively different setting for realistic atom-resonator 
couplings and loss rates, and is therefore observable in 
near future state of the art experiments involving just 
two coupled resonators. 

Conclusions - We have proposed the selective exci- 
tation of photonic many-body modes of interest in open 
resonator arrays using external driving lasers, over which 
we have full control of frequency and amplitude. We have 
shown how a combination of the equilibrium Hamilto- 
nian structure and non-equilibrium operation lead to an 
interaction-induced region of bunched emission. This fea- 
ture was found to persist in mesoscopic-sized arrays, and 
is also found under a more realistic array description, 
making its observation feasible in coming experiments. 
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